;+
; NAME:
;  hod_mcmc
;
; PURPOSE:
;  Generates MCMC samples for model parameters
;  using basic MCMC Metropolis-Hastings algorithm
; KEYWORD INPUTS:
;  nburn  - number of burn in samples (default 100)
;  nsamp  - number of MCMC samples for model parameters (default
;           1000)
;  nstep  - number of random walk steps to take when exploring model
;           parameter values (default 20)
;  savfile - set this to a file name that will save the parameters 
;            after every iteration
; KEYWORD OUTPUTS:
;  mean - mean values of parameters
;  error - 1 sigma errors on parameters
;  samps  - all the MCMC samples for each parameter (npars x nsamp
;              array)
;  lik - log likelihood values for each sample
; KEYWORD SWITCHES:
;  time - set to time MCMC
;  verbose - prints out various information as it progresses
;
;-
@hod_common
; The loglikelihood function
pro loglike,cube,ndim,npars,lnew
  COMMON hod_data
  COMMON hod_model
  
  ; first take the cube and convert from a uniform deviate 
  ; to a physical parameter value, drawn from the prior
  pars=cube
  nd=[ndim]
  for i=0,nd[0]-1 do begin
      pars[i]=hod_priorpar(cube[i],hod_priortype[i],hod_priorspec[*,i])
  endfor
   ; setting additional parameters to intial values
  if npars gt ndim then pars[ndim:npars-1]=initpars[ndim:npars-1]
  

  ; evaluate the model with these parameter values
  modval=hod_getmod(pars,hod_parfile)

  ; interpolate model to observed rp
  ; ** this isn't stritly needed now
  wp_mod=interpolate(modval.w_p_sum,findex(modval.r_p,rp_obs))
  
  ; evaluate the loglikelihood =chi2/2
  lnew=-total( ((wp_obs - wp_mod)/wperr_obs)^2 )/2

  if n_elements(dens_obs) gt 0 then if dens_obs gt -1 then lnew-= ((dens_obs - modval.density)/denserr_obs)^2 /2

  splog,'pars=',pars,' lik=',lnew
  ; store pars back in cube
  cube=pars
end
  


pro hod_mcmc,nburn=nburn,nsamp=nsamp,nstep=nstep,verbose=verbose,time=time,$
             mean=parmean,error=parerror,samps=samps,lik=liksamps,gibbs=gibbs,$
             savfile=savfile,resume=resume
   COMMON hod_model
   COMMON hod_data
 
   stime=systime(1)

  if keyword_set(resume) then begin
      splog,"Restoring previous MCMC runs from",savfile
      restore,savfile
      isamp0=isamp+1
  endif else begin
      isamp0=0

      npar=hod_npar
      ndim=hod_ndim
      
      if n_elements(nburn) eq 0 then nburn=100
      if n_elements(nsamp) eq 0 then nsamp=1000
      if n_elements(nstep) eq 0 then nstep=20 ; number of random walk steps
      
                                ; set up samps to hold the draws for sample
      samps=fltarr(npar,nsamp+nburn)
      usamps=fltarr(npar,nsamp+nburn)
      liksamps=fltarr(nsamp+nburn)
      step=fltarr(ndim)+0.5

      print,"Running MCMC.."
                                ; starting cube value chosen at random, calculate likelihood
      upars=randomu(seed,ndim)
      cube=upars
      
      loglike,cube,ndim,npar,loglik
      pars=cube
                                ; Now beginning the MCMC loop
      accept=intarr(ndim)
      reject=intarr(ndim)
      
  endelse


  steptime=systime(1)
  for isamp=isamp0,nburn+nsamp-1 do begin
     ; Now draw values for the model parameters using a random
     ; walk
      idim=0
      repeat begin
          for j=0,nstep-1 do begin
              ; storing current state
              prev_loglik=loglik
              prev_pars=pars
              prev_upars=upars
           
              ; step parameter within priors
              if keyword_set(gibbs) then upars[idim]=upars[idim]+randomn(seed)*step[idim] $
              else upars=upars+randomn(seed,ndim)*step
              upars=upars+(upars lt 0)-(upars gt 1) ; wraparound

              ; calculate likelihood for these parameters
              cube=upars
              loglike,cube,ndim,npar,loglik
              pars=cube
              
              ; Now decide whether to accept or reject this step
              dloglik=loglik-prev_loglik
              if dloglik gt 0 then begin
                  ; accept the step
                  accept[idim]++
                  if keyword_set(verbose) then print,'ACCEPTa!'
              endif else begin
                  logalpha=alog(randomu(seed))
                  if logalpha lt dloglik then begin
                      ; accept
                      accept[idim]++
                      if keyword_set(verbose) then print,'ACCEPTb!',dloglik,logalpha
                  endif else begin
                      ; reject
                      reject[idim]++
                      pars=prev_pars
                      upars=prev_upars
                      loglik=prev_loglik
                      if keyword_set(verbose) then print,'REJECT!',dloglik,logalpha
                   endelse
              endelse

           ; adjust the step size
              if keyword_set(gibbs) then begin
                 if accept[idim] gt reject[idim] then step[idim]=step[idim]*exp(1.0/accept[idim]) < 0.9999
                 if accept[idim] lt reject[idim] then step[idim]=step[idim]/exp(1.0/reject[idim]) 
              endif else begin
                 if accept[idim] gt reject[idim] then step=step*exp(1.0/accept[idim]) < 0.9999
                 if accept[idim] lt reject[idim] then step=step/exp(1.0/reject[idim]) 
              endelse
              
           endfor          
          idim++
      endrep until idim eq ndim or not keyword_set(gibbs)

      if keyword_set(verbose) then $
         print,"Draw ",isamp," pars=",pars," loglik=",loglik
      
      if keyword_set(time) then print,"draw ",isamp," took ",systime(1)-steptime," secs"
      steptime=systime(1)
      print,'****************************************************************************** '
      print,' '
     
      ;if isamp ge nburn then begin
          samps[*,isamp]=pars
          usamps[*,isamp]=upars
          liksamps[isamp]=loglik
      ;endif

      ; adjusting the stepsize
      if isamp gt nburn*2 then begin
         u=usamps[*,nburn:isamp]
         n=isamp-nburn+1
         m=total(u,2)/n
         step=sqrt( total( (u-rebin(m,npar,n))^2,2)/n )>1e-3
         ; resetting accept and reject rates
          accept=intarr(ndim)
          reject=intarr(ndim)
      endif
      splog,'step = ',step

      ; saving progress
      if n_elements(savfile) gt 0 then $
        save,filename=savfile,/all
 
  endfor

if keyword_set(time) then print,"Completed in",systime(1)-stime,"seconds"

; Now calculate summary statistics (and set model to mean)
parmean=total(samps[*,nburn:nsamp+nburn-1],2)/nsamp
parerror=sqrt( total( (samps[*,nburn:nsamp+nburn-1]-rebin(parmean,npar,nsamp))^2,2)/nsamp )
print,"Mean model parameters"
print,parmean
print,'1 sigma errors'
print,parerror    
        
        
  

end


  
  
